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Abstract 

The paper is devoted to penalty Robin-Robin domain decomposition 
methods (DDMs), proposed by us for the solution of unilateral multibody 
contact problems of elasticity. These DDMs are based on the penalty 
method for variational inequalities and some stationary and nonstation- 
ary iterative methods for nonlinear variational equations. The main result 
of the paper is that we give mathematical justification of proposed DDMs 
and prove theorems about their convergence. We also provide numeri- 
cal investigations of the efficiency of these methods using finite element 
approximations. 
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1 Introduction 

Contact problems of elasticity are widely used in many fields of science and 
engineering, especially in machine science, structural mechanics, geology and 
biomechanics. The brief overview of existing numerical and analytical methods 
for the solution of contact problems can be found in [TJ [2] . 

Efficient approach to the solution of multibody contact problems is the use 
of domain decomposition methods (DDMs). 

DDMs are well developed for the solution of linear boundary value problems, 
particularly for linear Poisson and linear elasticity problems El |5] . 

The construction of DDMs for unilateral contact problems, which are non- 
linear, are much more complicated. Among domain decomposition methods for 
unilateral two-body contact problems obtained on continuous level, one should 
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mention Signorini-Neumann [SJ [7], Signorini-Dirichlet [HI IS] and Signorini- 
Signorini [TU] iterative algorithms. All of these methods in each iteration require 
to solve a nonlinear one-sided contact problem with a rigid surface (Signorini 
problem) for one of the bodies, and a linear elasticity problem with Neumann 
[6j [7] or Dirichlet [U [9] boundary conditions on possible contact area for other 
body, or require to solve nonlinear Signorini problems for both of the bodies |10) . 
Moreover, to increase the convergence rate of Signorini-Dirichlet and Signorini- 
Signorini algorithms, it is recommended to perform an additional iteration, in 
which we have to solve linear elasticity problems with Neumann boundary con- 
ditions for both of the bodies [8j [10] . 

Domain decomposition method presented in work for two-body unilat- 
eral contact problem, is also obtained on continuous level. It is based on the use 
of augmented Lagrangian variational formulation and Uzawa block relaxation 
method. This domain decomposition method in each iteration require to solve 
linear elasticity problems for both of the bodies. 

On contrary, DDMs can be constructed on discrete level, after a discretiza- 
tion of corresponding continuous boundary-value problem. Among discrete 
DDMs for unilateral contact problems, one should mark out substructuring 
and FETI methods [H EH Efl E3 ■ 

In works [THl [TTJ HHJ HH] we have proposed on continuous level a class of 
penalty parallel Robin-Robin type domain decomposition methods for the solu- 
tion of unilateral multibody contact problems of elasticity. These methods are 
based on penalty method for variational inequalities and some stationary and 
nonstationary iterative methods for nonlinear variational equations. In each 
iteration of proposed DDMs we have to solve in parallel some linear variational 
equations in subdomains, which correspond to elasticity problems with Robin 
boundary conditions, prescribed on some subareas of possible contact zones. 
These DDMs do not require the solution of nonlinear one-sided contact prob- 
lems in each step. 

The main result of this paper is that we have proved theorems about conver- 
gence of proposed penalty Robin-Robin domain decomposition methods. The 
paper is organized as follows. In section 2 the classical formulation of multi- 
body contact problem in the form of the system of second order elliptic partial 
differential equations with inequality and equality constrains is given. In sec- 
tion 3 we consider the variational formulation of this problem in the form of 
convex minimization problem and in the form of elliptic variational inequality 
at closed convex set, and formulate theorem about a unique existence of the 
solution of this inequality. In section 4 we use the penalty method to reduce 
the variational inequality to an unconstrained minimization problem, which is 
equivalent to a nonlinear variational equation in the whole space. Later, we 
prove a theorem about a unique existence of the solution of the penalty vari- 
ational equation and a theorem about strong convergence of this solution to 
the solution of initial variational inequality. In section 5 we consider station- 
ary and nonstationary iterative methods for the solution of abstract nonlinear 
variational equations in reflexive Banach spaces. We prove theorems about con- 
vergence of these methods, and show that the convergence rate of stationary 
methods in some energy norm is linear. We also prove theorem about stability 
of stationary iterative methods to errors which may occur in each iteration. In 
section 6 we present parallel stationary and nonstationary penalty Robin-Robin 
domain decomposition methods for the solution of nonlinear penalty variational 
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equations of multibody unilateral contact problems. We prove theorem about 
convergence of these methods, and show that the convergence rate of stationary 
Robin-Robin methods in some energy norm is linear. In section 7 we provide 
numerical investigations of proposed domain decomposition methods using finite 
element approximations. The penalty parameter and mesh refinement influence 
on the numerical solution, as well as the dependence of the convergence rate 
of domain decomposition methods on iterative parameters are investigated. In 
conclusion section we summarize all results presented in the paper. 



2 Formulation of unilateral multibody contact 
problem 

Introduce the Cartesian coordinate system OX1X2X3 with basis vectors ei, e% 
and consider the problem of frictionless unilateral contact of N elastic bodies 
Q a C K 3 with piecewise smooth boundaries T a — dil a , a — 1,2,...,N (fig. I). 
Denote f2 = Ua=i 



Pa 




Fig 1. Unilateral contact of several bodies 

The stress-strain state in point x = (x\, x%, 23) T of each solid f2 Q is descried 
by the displacement vector u a (x) = u a , (x) e j , the symmetric tensor of strains 
Gj , and the tensor of stresses a a = a a y ej . These quantities 
satisfy Cauchy relations, Hook's Law and the equilibrium equations: 

3 

a ai j(x)= \] Cgijki(x)e a ki(x), x e Cl a , z,j = l,2,3, (2) 
fc,i=i 

E ^fo +Ui(*) = 0, xefi Q , i = 1,2,3, (3) 
where f a i (x) are the components of the volume forces vector f a (x) = f a , (x) e.; . 



3 



The elastic coefficients C a ijki(x) are measurable, symmetric, and bounded 
with constants 0<6<d<ooin the following sense [20] : 

3 3 3 

b M - H C ™i u M Sm i ( x ) £akl w - d 

(4) 

i,j=l ij.k.l—l k,l—l 

Suppose that the boundary T a of each solid consists of three parts: r™, 

r* s a , such that r Q = UrsU^, r^f]r-f]s a = 0, ^ 0, = rj», 

/Sq, 7^ 0. Boundary S a — U^eB Q S a p is the possible contact area of body fl a 
with other bodies, S a p is the possible contact area of body H a with body Hp, 
and B a C {1,2,..., AT} is the set of the indices of all bodies which are in contact 
with body H a . 

On each boundary T a let us introduce the local orthonormal coordinate 
system £ Q , £ a , n Q , where n a is outer unit normal to T a . Then the vectors 
of displacements and stresses on the boundary can be written in the following 
way: 

u„(x) = it Q £(x)£ Q + u qC (x)C q + w Q n(x) n Q , x e T a , 
tr a (x) = <r a (x) • n Q = CT a? (x)^ Q + cr aC (x)C Q + cr Qn (x)n Q , x e T Q . 

We assume that the surfaces S a p C T Q and Sp a C are sufficiently close 
[21] . Therefore n a (x) « — ng(x'), where x' = P(x) G is the projection of 
point x 6 5 Q( g on the surface S^q,. We denote by d a p(x) = ± ||x — x'|| 2 the 
distance in E 3 between bodies H a and Hp before the deformation. The sign of 
d a p (x) depends on a statement of the problem. 

On the part the kinematical (Dirichlet) boundary conditions are pre- 
scribed: 

u a (x) = z Q (x), x e T u a , (5) 
and on the part we consider the static (Neumann) boundary conditions 

<T a (x)=p a (x),xeC (6) 

where z a = z qC (x)£ q +z qC (x)C q + z q „(x) n Q and p Q = p a z(x)£ a +p a ((x.)C a + 
Pctn(x)n a are given boundary displacements and stresses. 

Further, for the simplicity of variational formulations and proofs, we assume 
that all bodies are rigidly fixed on the surface rjj, i.e. 

z a (x) = 0, xeC (7) 

On the possible contact areas S a p, a = 1, 2, ...,N, (3 £ B a the following unilat- 
eral contact conditions hold: 
absence of extension 

<7 an (x) = cr (3 „(x') < 0, (8) 

absence of friction 

<t q5 (x) = crpe(x!) = 0, cr Q<; (x) = cr^ c (x') = 0, (9) 
mutual nonpenetration of the bodies 

u an {x) + up n (x') < d a p(x), (10) 
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contact alternative 

(it a „(x) + u i s n (x') - d Q(3 (x) ) cr Q „(x) = , (11) 

where x G S^, x' = P(x) € S' ( g Q ,. 

The system of second order partial differential equations (HJ) - © with 
boundary conditions ([5]) - (1111) is the mathematical formulation of frictionless 
unilateral multibody contact problem of elasticity. 

Note that the contact problem (jXJ) -■ fl3J), © - (|TTj) is nonlinear, since the 
real contact areas are unknown. 



3 Variational formulation of the contact prob- 
lem 

Let us give the weak formulation of the contact problem ([I]) - © , © ~ HQ hi 
the form of variational equality and convex minimization problem. 

For each body £l a , a = 1,2, ...,7V consider Sobolev space V a = [i? 1 (17 Q )] 3 

with scalar product (u a , ^ a ) Va = Yn=i Ia a (u a iV a i + Y?j=i T^T^f) ^> 
€ V a and norm HuqH^ = J (u Q ,u a )^, u Q G V a . 
Introduce the following closed subspace in V a ; 

V° = {u a : u^K, lrK)=0onP o }, (12) 

where Tr^ : V a — > [i/ 1 / 2 (r^)] 3 is surjective, linear and continuous trace op- 
erator [201 HI] • Space V® is a Hilbert space with the same scalar product and 
norm as in V a . 

Consider the space Vo, which is the direct product of spaces V£: 
Vo = V? x ... x y° = {u = ( Ul , u N f , u a e V°, a = 1, 2, AT I , (13) 

and define in it the scalar product and norm: (u,v) Vo = Y] a —i (u Q ,v Q )y , 
IjuH^ = ^/(u, u)^, u, v € Vo- Note, that Hilbert space Vq is closed reflexive 
Banach space. 

Now, let us introduce the closed convex set of all displacement vectors in Vq 
which satisfy nonpenetration contact conditions (1101) : 

A' = {u: ueV , u an + up n < d aj3 on S a p , {a, P} £ Q} , (14) 

where Q = { {a, f3} : a € {1, 2, A^} , /3 G B a } is the set of all possible un- 
ordered pairs of the subscripts of bodies which are in contact with each other, 
and d a0 e ffo /2 (£ a ), {a, /?} € Q, E Q = int (r a \r£), a = 1, 2, N. 

The quantities u Q „, a = 1, 2, JV in (1141) have to be understood in the 
following way 

u an = n a ■ Tr°(u Q ), u Q e V„, 

where Tr° : V® — > [Hqq 2 (T, a )] 3 is surjective, linear and continuous trace oper- 
ator on surface E Q = int (T a \r^) HUES], and n Q € [i 2 (S Q )] 3 . 

Note, that all equalities and inequalities in spaces L2, if 1 / 2 , H^ 2 and H 1 
hold almost everywhere. 
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Since set K is a closed convex subset of Hilbert space Vo [101 HD 122] , it is 
weakly closed [23] . 

In space Vo x Vo consider bilinear form A(u,v), which represents the total 
deformation energy of the system of bodies: 



A(u,v) = ^2 a Q (u Q ,v Q ), u,veVo, (15) 

a=l 

«c t (u a ,v )= / a a (u a ) : e a (y a ) dtt, u a ,v a eV°. (16) 
In space Vo define linear form L (v) , which is equal to the external forces work: 



N 



L(v) = Y,la(Va), v€Vq, (17) 



a=l 



Uv a )= [ f a -v a dQ+ [ p a - Tr° a (w a )dS, v Q GV Q . (18) 

where f a G [L 2 (0>)] 3 , p Q G [ffoo 1/2 ( s «)] 3 , a = 1,2,..., TV. 

Lemma 1. If T a = dft a , a = 1,2,...,N are piecewise smooth, ^ 0, 

r« = f„ G [L2(fi Q )] 3 , Pa G [-^oo l72 (^a)] 3 ' an< ^ condition holds, then the 
bilinear form A is symmetric, continuous and coercive, and the linear form L 
is continuous, i.e. 

(Vu,vGVo){A(u,v)=A(v,u)}, (19) 
(3M > 0) (Vu, v G Vo) { |A(u,v)| < M \\u\\ Vo \\v\\ Vo } , (20) 
(3B > 0) (Vu G Vo) {A (u, u) > B \\u\\ 2 Vo } , (21) 

(3T>0)(VvG V ){\L(v)\ <T\\v\\ Vo }. (22) 

Proof. The symmetry of bilinear form A (u, v) follows from the symmetry of 
elastic constants in the Hook's law ©, and its continuity and coercivity follows 
from the continuity and coercivity of bilinear forms a a (u Q , v Q ) for each body. 
The continuity of L (v) follows from the continuity of linear forms I a (v a ), which 
can be proved using trace theorems [3U] . The detailed proof of this Lemma can 
be found in [T5]. □ 

According to [2T] [23], the original contact problem © - ©, © - $TQ 
has an alternative weak formulation as the convex minimization problem of the 
quadratic functional on the set K: 

F(u) = ±A(u,u)-£(u)^nun. (23) 
z u ek 

Using the general theory of variational inequalities [HJ HI] we have proved 
in [21] [23] the next theorem. 

Theorem 1. Suppose that bilinear form A and linear form L satisfy prop- 
erties Aiyp - i22f) , and set K is convex closed subset of Hilbert space Vq . Then 
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the minimization problem \23\) has a unique solution u G K , and this problem 
is equivalent to the following variational inequality: 

F'(u, v - u) = A (u, v - u) - L (v - u) > 0, VveiC. (24) 



4 Penalty variational formulation of the prob- 
lem 

To obtain a minimization problem in original space Vb, we apply a penalty 
method 24, 25 to convex minimization problem (1231) . 

For the violation of nonpenetration conditions we use the penalty in the 
following form |T71 [20] : 

Je(u) = — ^ / {d a p - u an - Uf3 n y dS, (25) 

where 9 > is a penalty parameter, y~ = min{0, y}. 

Let us consider the following minimization problem with penalty in space 

V : 

Fg{u) = -A (u, u) - L(u) + J fl (u) -> min . (26) 

Z u £Vb 

Note, that the introduction of penalty corresponds to the introduction of condi- 
tional intermediate Winkler layer between bodies with stiffness coefficient 1/6. 
The quantity o- a p n = cF an — n = (d a p — u an — up n )~ /0 has a sense of the 
normal contact stress between bodies fl a and Hp, and the penalty Jg(u) repre- 
sents the total normal contact stresses work. 

Now consider the properties of penalty term (|25p in more detail. Functional 
Je(u) is nonnegative 

(Vue V o ){J e (u)>0}, (27) 
and Gateaux differentiable in Vq: 

J' (u,v)=~ ^2 ( {d a p - u an - Uf3 n )~ [v a n+vp n )dS, (28) 
{ a , /3}eQ Js "f 

Furthermore, the Gateaux differential J' e (u, v) is linear in v. 

Lemma 2. If surfaces S a p, {a, /?} G Q are piecewise smooth, and d a p G 
l li 

Hqq (S q ), then Jg(u, v) satisfy the following properties: 

(Vu G Vb) (aiZ > 0) (Vv G Vo) { |4(u,v)| < R\\v\\ Vo } , (29) 

(Vu, v G Vb) { J' e (u + v, v) - J' e (u, v) > 0} , (30) 

(3D > 0) (Vu,v,w G Vo) { |4(u + w,v) - 4(u,v)| <D\\w\\ Vo ||w|| Vo } . (31) 

Proof. It is obvious to see, that functional J' g (u,v) is linear in v. 
At first, let us show the satisfaction of property (l29l) . Let us write J#(u, v) 
in the extended form: Jg(u,v) = Y,{ a , ^ e Q j'apfa v ), where 

iU u ' v ) = 4[/ s g a p(u) n a -Tr° a (v a )dS + g a p{u) ■ Tr^(v^) dsj , 

(32) 
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and g a p(u) = {d a p - u an - up n ) . 

Taking into account the following inequality for real numbers: 



5Z C M ^ m ^2 c t> c *eR, i = l,2,...,m, men, (33) 

\i=l / i=l 

and Schwarz inequality, we obtain 

^jf g a p(u) n a -Tr° a (v a )ds\ < 3 ?a/} (u) \\ Tr°( v a )|| [ i2(Sa/3)]m , (34) 



where g aj g(u) = ||g Q /j(u) n a || [i2(s ^ 3 + e Q /3, e a( g > 0, /3 G B a , a = 1,2,...,N. 
From trace theorems in [2D], it follows, that 

{3T aP > 0) (Vv Q G V°) { ||Tr°(v a )||j L2(s ^ )]3 < 1% Wf Va ) ■ (35) 

Substituting (l35|) into (|34|) . we come to inequality / s ^ g a p{u) n a ■ Tr° ( v Q ) < 

s Q/3 (u) Hvally^, where s a0 (u) = T a/3 y / 3g^g(u) > 0, /3 G S Q , a = 1,2, AT. 

Similarly to this, we obtain the same inequality for the second term of rela- 
tionship (|3"2"|) . Hence 

|ja/}(u,v)| < ~ (s Q/J (u) ||v a || Va +S/3a(u) 1 1 V,3 1 1 V;j ) • 

As a result, we find 

I4(u,v)|< K^u,v)|<i?(u)||v|| yo , 

{a, 13} £Q 

where .R(u) = | X)r Q| /3} e q (*a^( u ) + «/3a(u)) > 0. Inequality fl2S) is proved. 
Now, we prove that condition (|3"0|) holds. For this we use the next inequality 

(Vy,zGE){ [(y-z)--y-] z < 0} , (36) 

Rewrite Jg(u + w, v) — Jg(u, v) in the following way: 

4(u + w,v) - Je(u,v) = -- /ia/3(u,w,v), (37) 

{a, /3}eQ 

where /i Qj g(u,w, v) = J s r Qj g(u, w) ( ) dS, r a/3 (u,w) = g a ^(u+w)- 

5a/3(u), x G S a p. 

In view of property (|36p . we obtain 



(Vu, v G Vb) {r Q/3 (u, v) (u an + vp n ) < 0, x G S Qj a} . 
Therefore, since > 0, we come to inequality 

J' g (u+v,v)-J' g (u, v) = -~ / r a/3 (u, v) (v Q „ + u^n) dS > 0, Vu, v G V . 
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Finally, let us prove the satisfaction of property (f5T|) . Using the next in- 
equality for real numbers 



we obtain 



(Vy,zeR){\y- -z~\ < \y - z\] 



r Q/3 (u, w) < \w an + Wp n \ , x G S a p. 



(38) 



(39) 



Taking into account, that the components of outer unit normal to <9f2 Q satisfy 
property max \n a j \ < 1, and using inequalities (|33p and (f3T)|). we hnd 

/ r^(u,w)dS < 6 (||Tr Q (w Q )||^ ( ^ )]3 + ||Tr°(w^)||^ ) . (40) 

Now, let us write h a /3(u, w, v) as follows: 

h a p(u, w, v) = / r a (s(u,w)v an dS + / r aj3 (u, w) i>g „ dS. 

Consider the term J s ^ r a ^(u, w) v an dS in more detail. Let us use inequalities 
(|3"3"| . (HO]), and Schwarz inequality: 

/ r a fs(u 7 w)v an dS \ < / ^(UjWjdyS / v 2 an dS < 



<18(||Tr°(w Q )||^ ( ^ )]3 + ||Tr«(w,) 
Using inequality (|35[) . we find further 



K(V Q )| 



J r a/3 (u,vf)v an dS < T a/ 3i IjVally^ Q|w„|| Voi + llw^ll^^ 



where T a p\ = 2>T a pT*o^/2 > 0, T*^ — max{T a p, Tp a } > 0. In much the same 
way we come to inequality 



/ r a p(u,w)vp n dS 

JS„R 



< T al 32 \\vp\\ V/j [ \\yr a \\ Va + ||w^|| 



where T Q/32 = 3 T$ a T*^ > 0. 

Taking into account the last two inequalities, we establish 



\h a p(u,w,\)\ < 



/ r Q/3 (u,w) v an dS + / r a j 3 (u,w)v l3 



dS 



< 



< 



C al 3 (\\v a \\ Va \\vf a \\ Va + \\Va\\ Va \\yrp\\y e + \\vp\\ Vft \\Wa\\ Va + \\^\\v ff ll W ^llv^) 



where C a p = max {T Q/3 i, T Q/32 } > 0. 
As a result, we obtain 



1 N 

\J'e(u + w,v)-J' g (u,v)\<- Kp{u,w,v)\ < 



<*,P=1 
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2C I N N N 

a=l a=l /3=1 



^ o - Ell^ilv.ElKHv^-PlMI HI, Vu,v,weV„, 



o=l 0=1 

where D = 2C (N + 1) /(9 > 0, C = max C* Q/3 > 0. □ 

l<a,P<N 

Theorem 2. Suppose that Vo is closed reflexive banach space, Jg(u) is 
Gateaux differentiable, and conditions fWj) - HUtty . [Wfy , fj§|) - K31)) hold. Then 
there exist a unique solution of nonquadratic minimization problem \26\) in Vo, 
and this problem is equivalent to the following nonlinear variational equation: 

iftu.v) = A(u,v) + 4(u,v)- J L(v) =0, VveVb, ueV . (41) 

Proof. As shown in [T7], due to the properties (TH))) - (j2"2"|) . the functional 
F(u) is strictly convex, and coercive ( lim F(u) — 00), and the differential 

F'(u, v) is linear and continuous in v. 

From the property (I5P1) . it follows that the penalty term Je(u) is convex in 

Now consider the properties of functional 

Fg{u)=F{u) + Jg{u), ueV . 

This functional is Gateaux differentiable in Vo: 

F e '(u,v) = F'(u,v)+4(u,v), u,veVo, 

and strictly convex, as the sum of convex functional Jg(u) and strictly convex 
functional F(u). In addition, since functionals F'(u,v) and Jg(u, v) are linear 
and continuous in v, it follows that Fg(u, v) is also linear and continuous in v. 
Hence, according to [24], functional Fg(u) is weakly lower semicontinuous. Due 

to the coercivity of F(u) and property (26), we obtain that lim Fg(u) = 00. 

Il u ll v- — 

Since Vb is closed reflexive Banach space, Fg (u) is weakly lower semicontin- 
uous and coercive, then according to the theorem of existence of the minimum 
in reflexive Banach space |24) , there exists a solution of minimization problem 
(l26l) in the space Vq. 

Furthermore, the functional Fg(u) is strictly convex and Gateaux differen- 
tiable in Vq. Therefore, according to the theorem about necessary and sufficient 
conditions of the minimum in reflexive Banach space |24j . the solution of prob- 
lem (1261) is unique, and this problem is equivalent to the variational equation 

dm. □ 

Now let us prove that the solution of penalty variational equation (|41l) con- 
verges strongly to the solution of original variational inequality as 9 — > 0. 
Let us rewrite variational equation (|4ip in the following equivalent form 

fi(u,v)=i(u,v)-(i,v) + l(*(u),v)=0, VveVo, ue Vo, (42) 
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where (Y, u) = Y(u) is the action of a functional Y £ Vq* on an element u £ Vb, 
VJj* is the space dual to Vq, $ = ^' : Vb — > Vq* is Gateaux derivative of functional 
*(u) = 0Jg(u), and (*(u),v) = (*'(u),v) = 9J' e (u,v). 
We have proved the next lemma. 

Lemma 3. Suppose that properties L30\) and L31\) hold. Then operator $ : 
Vq —} V * in problem is a penalty operator for the kinematically allowable 
displacements set K , i.e. 

1) . $ is monotone in Vq: 

(Vu,ve V ){($(u)-$(v),u-v) >0}; 

2) . $ satisfies the Lipschitz condition in Vq: 

(3C > 0) (Vu,vg V ) { ||$(u) - $(v)|| Vo . < <7||u - vllvj ; 

3) . The kernel of operator $ is equal to set K : 

Ker ($) = {u : u G Vq, $(u) = 0} = K. 

Proof. The monotonity of operator $ follows from the condition (|30|) , and 
the satisfaction of Lipschitz condition follows from the property (f3~Tj) . 

If u € K, then $(u) = 0, and, on the contrary, if $(u) = 0, we have u G if. 
Hence, Ker($) = if. For more details see [T71 [T^] - □ 

Now, using the results of works [251 121)1 \Tf\ . let us prove the proposition 
about the convergence of penalty method, applied to the variational inequality 

HMD- 
Theorem 3. Suppose that K C Vb is a convex closed subset of Hilbert 
space Vb, bilinear form A (u, v), u, v G Vb and linear form L (v), v G Vb satisfy 
properties h!9\) - 1221) . : Vq —¥ V * is penalty operator for set K, u £ K 
is a unique solution of variational inequality \2J$ , and Ug £ Vq is a unique 
solution of penalty variational equation \4ty with penalty parameter 9 > 0. Then 
ug — > u stronqly in Vq, i.e. \\vlq — fill,, — > 0. 

Proof. In works [IS [53] it is proved that, if conditions (JH| - © hold, $ 
is a penalty operator for set K and there exist solutions of problems (|24|) and 
(|42l) . then the sequence {u#} is bounded: 

(a(7G (0; co)) (V0>O){||u e || yn <c), 

and there exists such subsequence {u^} C {u#}, which converges weakly in Vq 
to some solution of variational inequality (|24[) . i.e. 

(3{u Sl } c {u e }) (VY G v*) | (y,a fll ) ei ^ (^Q>} ■ 

Moreover, in [25, 26. it is shown that any weakly convergent subsequence of 
sequence {ug} converges weakly in Vq to some solution of variational problem 

Now let us assume that variational problems (jM)) and (|42l) have unique 
solutions. Then, as follows from above, the sequence {ug} has a unique partial 
weak limit u G K. 
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Since the sequence {ug} has a unique weak limit point, and is bounded, 
then according to theorem in [24] . it is weakly convergent to this point, i.e. 



where u £ K is the unique solution of variational inequality (f2"4"|) . 

Further, let us show that {u«} converges strongly to u £ K as 6 — » 0. 
Due to (g3]), we get 

(L,u e -u) -> 0, A(v,u e -u) = K(v),u e -u) -> 0, VveV , (44) 

where ^4'i( v ) is the Gateaux derivative of functional A±(v) — \A (v, v), v £ Vo- 
Since u $ is a solution of penalty variational equation (|42l) , it is obvious that 

(L,u 9 - v) = A(u g ,u e - v) + i ($(u 8 ),u e - v) , Vv G if. 

f 

Taking into account the monotonity of penalty operator <£> and the property 
(Vv G K) {$(v) = 0}, we obtain 

(L,U0 - v) = A (u e , u e -v)+i ($(u e ) - $(v),u 9 - v) > ^4 (u 8 , u e -v), Vv G if. 

In view of this property and the nonnegativity of bilinear form A, we get the 
inequality A (u, flj- u) — (L, u 8 — u) < A (tig, u.e — u)~(L, xig — u) < 0. Hence 



Passing to the limit in expression (|4"5]l as — >• 0, and taking into account 
property (|44|) . we obtain 



Further, in view of coercivity of bilinear form A, it follows that 
< B ||u* - u||y < i4 (u fl , u - u) — A(u,u e - u) -> 0, B > 0. 

As a result, we establish that IIua — ull,, — > 0. □ 

Thus, using the penalty method we have reduced the solution of original 
variational inequality (|24p to the solution of the nonlinear variational equation 
in the whole space Vo, which depends on the penalty parameter 8 > 0. We also 
have proved the existence of the unique solution of penalty variational equation 
(|4"2"]l and its strong convergence to the solution of initial variational inequality 
(|24| as penalty parameter tends to zero. 

5 Iterative methods for nonlinear variation equa- 
tions 

Consider an abstract nonlinear in u variational equation in form (|4"Tj) . where 
Vq is closed reflexive Banach space, A (u, v) is bilinear form, set in Vq x Vo, 




(43) 



A(u,Ug - u) < A(ug,Ug - u) < (L,Ug - u) . 



(45) 



A(ug,Ug - U) ->• 0. 
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L (v) is linear functional, and the term J' s (u, v) is linear in v and nonlinear in 
u. Suppose that conditions - (j2"2")l . (|2T))) - (f3"Tj) are satisfied. Hence, there 
exists a unique solution of the problem (jlTj) . 

For the numerical solution of the nonlinear variational equation (|4ip let us 
use the following iterative method [TBI H71 fTHl fT9] : 

G (u fc+1 , v) = G (u fc , v) - 7 [A (u fc , v) + J> fc , v) - L (v)] , k = 0, 1, ... , (46) 

where G(u,v) is some bilinear form assigned in Vq x Vq, u k £ Vq, k — 1,2,... 
is the fc-th approximation to the exact solution u £ Vq of the problem (l4l"|) . 
u° e Vb is an initial guess, and 7 € M is an iterative parameter. 

We have proved the following proposition about the convergence of the iter- 
ative method (|4"6")l . 

Theorem 4. Suppose that bilinear form G (u, v) is symmetric, continuous 
and coercive: 

(Vu,ve V ){G(u,v) = G(v,u)}, (47) 
(lM > 0) (Vu, v £ V ) { |G(u,v)| < M||u|| Vo \\v\\ Vo } , (48) 

(BB > 0) (Vu £ V ) {g(u,u) > SHull 2 ^} , (49) 

properties H2U\) - h2'2\) , \29fl - \31}) are satisfied, and the iteration parameter lies 

in the interval 7 £ (0; 72), 72 = 2BB j Ml , M» = M + D. 

Then the sequence {u fe }, obtained by iterative method J^6] ) converges strongly 
in Vq to the exact solution u £ Vq of the variational equation J^J[ ), 
i.e. ||u fc — u||y o — > 0, anii t/ie convergence rate in the energy norm 

|ju|| G = y/G (u, u) is linear: 



||u fc+1 -u|| G <g||u fc -u|| G , q = Jl- 1 (2B- 1 M%/B) /m < 1. (50) 

Moreover, the maximal convergence rate reaches as 7 = 7 = -BB . 

Proof. Since bilinear form G (u, v) is symmetric, continuous and coercive, 
we may introduce a scalar product and norm 

(u, v) G = G(u,v), ||u|| G = ^G(u,u), u,v£V . 

From properties (|48|) and (|49l) . it follows that norms || • || G and || • ||y o are 
equivalent in space Vq. 

In each step k £ {0, 1, ...} of method (|4"6")l . we have to solve the linear varia- 
tional problem: 

G(u,v) = F fc (v), VveVo, ueV , (51) 

where Y k (v) = G (u k , v) — 7 [A (u k , v) + J' g (u k , v) — L (v)l is linear in v, u fe e 
Vo- Using properties (HOI), ([221), O and (|i5j). we obtain that y fe (v) is contin- 
uous: 

(3z fc >o)(Vvey ){|r fc (v)| <z fc ||v|| yo }, (52) 

where Z fc = M H^H^ + | 7 | (m \\u k \\ Vo + R (u fc ) + Tj + e > 0, e > 0. 

Since conditions (|4T)l - (|4"9"|) and (1521 are satisfied, we see that the problem 
([5Tj) has a unique solution u = u fe+1 e Vq. 
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Now let us show, that the sequence of solutions of problems (|5ip converges 
strongly to the solution of initial variational equation (I41[) . 

Suppose that u 6 Vb is the exact solution of problem (pfTj) . Introduce a 
notation <^ fc := u fc — u € Vb, fc = 0, 1, and rewrite ([15} as follows: 

G (u + <p k+1 , v) = G (u + y> fc , v) - 7 [A (u + <p fe , v) + J' g {u + <p k , v) - L (v)] . 

Subtracting from this expression the identity G (u, v) = G (u, v)— -f[A (u, v)+ 
Jg(u, v) — L (v)], we obtain 

G (y> fe+1 , v) = G ( V fe , v) - 7 [A (<p k , v) + 4(u + p*, v) - 4(u, v)] . (53) 

Let us define a functional Hg(u, w, v) = Jg(u+w, v) — J' 8 (u, v), u, v, w g Vb, 
which is linear in v. 

Due to properties (j30|) and (|31j). the following conditions hold: 

(Vu,ve Vb){i/ e (u,v,v) >0}, (54) 

{3D > 0) (Vu,v,w e Vb) { |ff fl (u,w,v)| < D \\w\\ Vo \\v\\ Vo } . (55) 
Let us rewrite expression (|53[) in the form 

G(<p k+1 -<p k ,v) = - 1 [A(<p k ,v) + H e (u,<p k , v)], VveVb. (56) 

If we take v := ^ fc+1 — y> fe in ([56|h we shall have 

\\<p k+1 ~<p k \\ G < \l\ (\A(<p k ,<p k+1 -ip k )\ + \H e (u,<p k ,<p k+1 -tp k )\) . 

Taking into account the continuity of bilinear form ([%(!)) . and property (|55l) . 
we get 

\\(p k+1 -<p k \f G < |7|M» 11^*11^ ||^ fc+1 -v? i: || Vo , M* = M + £> >0. 
Further, in view of relation between norms 



\<P h+1 -V k \\ Vo < \\<p k+1 -<p k \\ G /VB, (57) 



we come to inequalities 



1 h k+1 - f k \\ Vo < \\<p k+1 ~ f\\ G < ^ h k \\ Vo ■ (58) 



Now let us take v := <p k+1 +tp k in expression (|56"|) . Then G (y> fe+1 — y> fc , (p k+1 + 
tp k ) = — 7 [A (ip k ,ip k+1 + (p k ) + Hg(u, (p k ,(p k+1 + ip k )] . This relation can be 
written as 

\W k f G - \W k+1 \\ 2 G = 7 + 2Hg(u,ip k ,ip k )+ 

+A (<p k ,<p k+1 - tp k ) + H e (u,ip k ,ip k+1 ~<p k )]. 
In view of properties (|20[) and (|55l) , we come to inequality 

A{<p k ,<p k+1 -<p k ) +H B (u,<p k ,<p k+1 -<p k ) > -M* \\<p k \\ Vo \\<p k+1 -<p k \\ Vo . 
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Suppose that 7 > 0. Then, taking into account the coercitivity of bilinear 
form ([2"T]) . property (I5~1T) . and the previous inequality, we obtain 



W k \\r-\W k+1 \\a>l 



2B h k \\ Vo ~ M AW k \\ Va h k+1 -* k \\v 



In view of inequalities (|57[) and (|58|l . we find further 

|k fc ||c-|k fe+1 ||c>^(2i3-7A4 2 /s) (59) 
If the following inequality 
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{^B-^Ml/B) >0 (60) 



holds, then the sequence ||y fc || G . will be monotonically nonincreasing: ||¥> fc || G > 
||y> fc+1 || 2 , and, hence, Hv^ll^ — > w, where to > 0. Passing to the limit as 



fc — > oo in expression (j5"9")l . we obtain > jfe y2B — jM 2 j Bj w, i.e. u> = 0, 
and, therefore — > 0. From inequality (|B(J1) we establish the interval of 

^ fc— s-oo 

allowable values of iteration parameter 7: 

7G (0; 72 ), 12 = 2Bb/mI. 
Since the norms || • || G and |j • || y are equivalent, we have ||y fe || 1/ — > 0, and, 

II II Vo fc—^oo 

hence, ||u fc — u|L, — > 0. 

Using inequality (|59[) . we find estimate 

|k fe+1 || G <g 2 ||^H G , g 2 (7) = l-^7+^7 2 - (61) 
11 llG 11 IIG u; M MB 

It is not hard to show, that q 2 £ (0; 1) for 7 G (0 ; 72). We obtain this from 
the next relations: 

q 2 {0) = g 2 ( 72 ) = 1, 7 = arg min g 2 ( 7 ) = 72 /2 = Bb/m 2 . 

76(0;72) ' 

As follows from estimate (|6ip , the convergence rate is maximal if the param- 
eter q is minimal, i.e., if 7 = 7. □ 

Remark 1. Suppose that term J'g{u, v) is Gateaux differentiable in u. Then 
conditions \29\) , in theorem 4 can be replaced by the following properties 

(Vu,veF ){4'(u,v,v) >0}, (62) 

(323 > 0)(Vu,v,w £ V ) { |4'(u,v,w)| < I>||v|| Vo ||w|| Vo }. (63) 
Proof. Let us apply to J'g'(u, v, w) the Lagrange formula of finite increments 

(Vu,v,w G Vb)(3r G (0;1)){ J^(u + rw,w,v) = J^(u + w,v) - 4(u,v)}. 

Then the satisfaction of conditions (|6"2")l and (|6"3"|) yields the satisfaction of prop- 
erties dnSJ and (EJ). □ 
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Now let us investigate the stability of iterative method (PfoT) to errors. 

Suppose that conditions of theorem 4 are satisfied. Then for any u fc 6 Vq, 
k = 0, 1, ... there exists a unique solution u = u fe+1 e Vb of problem ([ST]) . 

Therefore, there exists an operator R : u k G Vb — > u fc+1 e Vb, which maps 
every element u k E Vq to the solution u = u k+1 of the problem (|5ip . and the 
iterative method (|46l) can be written in the following form 

u fc+1 = R(u k ), k = 0,1,.... (64) 

Now assume that in each step k of iterative method (|6"4"]l we get some com- 
putational errors. Then this iterative method will take the form: 

u°=u n +e°, (65) 

u fc+1 =R{u k )+e k+1 , A: = 0,1,..., (66) 

where ii fc+1 , A: = 0,1,... is an approximate solution of problem (|5T]) . £ fc+1 , 
A; = 0, 1, ... is a computational error, which occurs in each step k, and e° is an 
initial guess error. 

Corollary 1. Suppose that conditions of theorem 4 are satisfied, and 
errors which occur in each step k of iterative method |^6] ) are uniformly bounded, 
i.e. 

(3e>0) (Vfc € {0, 1, ...» {\\e k \\ G <e}. 
Then the following estimates hold: 

||ii fc - u|| G < q ||u fc_1 - u|| G +£, (67) 
||ti fe -u|| G < 9 fe ||u -u|| G + ^-, (68) 
||ti fe ~ u|L < ||u fe - il^lL + -i-, (69) 

|[ti fc -u*|| G < 2^11^-01^ + ^, fc = l,2,..., (70) 

where u € Vb is the exact solution of problem 

The proof of this proposition follows from property (1501) . 

Thus, the errors, which occur in each step of iterative method (|46p . do not 

accumulate. 

Now consider the nonstationary iterative method for the solution of nonlinear 
variational equation (|4ip , where bilinear forms G (u, v) are different in each 
iteration (191158]. 

In space Vb x Vb introduce a sequence of bilinear forms {G k : Vb x Vb — > M} , 
k = 0, 1, which satisfy the property 

(V Y G Y *) (3 ! u e Vb) (Vv g V ) {G fe (u, v) - Y(v) = 0} . 

For the solution of nonlinear variational equation (|41l) , we have proposed the 
following nonstationary iterative method [TH] [5S] : 

G*( u *+i )V ) =G fc (u fc ,v)- 7 [A(u fc ,v) + 4(u fe ,v) -L(v)] , fc = 0,l,..., 

(71) 
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where u fe G Vq is the fc-th approximation to the exact solution of problem ([4"Tj) , 
and 7 £ R is an iterative parameter. 

We also have proved the next proposition about the convergence of this 
method. 

Theorem 5. Suppose that conditions \20j) - H22\) . h29\) - \31}) hold, bilinear 
forms G k (u, v) satisfy properties 

(Vu,v £ V ) |G fe (u,v) = G fe (v,u)} , (72) 

(BM > 0) (Vfc £ {0, 1,...}) (Vu, v € Vb) { |G fe (u,v)| < M||u|| % H^} , (73) 

(3B > 0) (Vfc £ {0, 1, ...}) (Vu £ V ) {G k (u, u)>B \\uf Vo ] , (74) 

(3fc G {0, 1,...}) (Vfc > fc ) (Vu G V ) {G k (u,u) > G k+1 (u,u)} , (75) 

and iterative parameter 7 lies in the interval 7 G (0; 72), 72 = 2BB j (M + D) 2 . 

Then the sequence {u k } obtained by the nonstationary iterative method (7l\ ) 
converges strongly in Vq to the exact solution u £ Vq of the variational equation 
OB , i.e. \\\i k - ull -> 0. 

The proof of this theorem is similar to the proof of theorem 4. We omit 
the details. 

Note, that in iterative methods and (fTTj) we can also take the parameter 
7 differently in each iteration. The purpose of the nonstationary choice of 7 
might be the improvement of the convergence rate of iterative methods (|46[) 
and (|7T|) . 



6 Parallel domain decomposition schemes 

Note, that in the most general case iterative methods (|4"6"]l and (fTTj) applied to 
solve nonlinear penalty variational equations (|41l) for multibody contact prob- 
lems do not lead to domain decomposition. Therefore, we now consider such 
variants of these methods, which lead to domain decomposition, namely, which 
reduce the solution of the original multibody contact problem in to the so- 
lution of a sequence of separate linear variational problems in subdomains tt a , 
a = 1,2,..., N. 

Let us take the bilinear form G in the iterative method (|4"6l as follows [T51I19] : 
G(u,v) = ^(u,v)+X(u,v), u,veVb, (76) 
where X (u, v) : Vq X Vq — > R is the next bilinear form [18l [19] : 

X(U,V) = - ^ J (UanVanlpa/3 + Up n V p n 1p p a ) dS, U,VGVb- (77) 

Here ip a p(x-) — |o, x G S' Q ^\S , ^ ( g| V |l, x G S^g | are characteristic functions 

of some given subareas S^p C S a /s of possible contact zones S a p, a = 1,2, N, 
p£B a . 
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The iterative method PS|) with bilinear form ([77)1 can be written in the 
following equivalent way: 

A(u k+1 ,v) + X(u k+ \v) =L(v)+X(u fe ,v)-4(u fe ,v), Vv e V Q , (78) 

u k+1 = 7 u fe+1 + (l- 7 )u fc , fc = Q,l,.... (79) 

Lemma 4. Suppose that surfaces S a p, {a,/?} € Q are piecewise smooth. 
Then bilinear form PTa) is symmetric, continuous and nonnegative, i. e. 



(Vu,vey„){X(u,v) = I(v,u)}, (80) 
(3Z > 0) (Vu, v e Vb) { |X(u,v)| < ZlHIvj ||v|| v . } , (81) 

(Vu € V ) {X (u, u) > 0} . (82) 

Proof. It is obvious, that conditions ([5U)) and (|82p hold. Thus, let us show 
the continuity of bilinear form (1771) . 

We can write X (u, v) = J2{ a , ^eQ^Kv), where 



dS , {a,£}eQ. 



X Q ^(u,v) = i^y Tpal3U an V an dS + J 

The first term can be written in the following form: 

/ 4>afiu an v an dS = / ip a (i [n Q ■ Tr°(u Q )] [n Q ■ Tr° ( v„)] dS, 

Taking into account, that the functions tp a p and the components of unit normals 
n a are bounded, and using inequality (|33[) and Schwarz inequality, we obtain 



dS^j <9||Tr° Q (u Q )| 
In view of inequality f|35[) . we find further 



2 \\Tr° fv ^lll 2 



ip a pu an v an dS < 3 W a/3 ||u a || v ||v a || v , W a p=T*g+e a p>Q, e a p > 0. 
Similar inequality can be obtained for the second term of X a p. Thus 

|X a/3 (u,v)| < (l|uc|lv a ll v allv a + \\ u p\\v p H^HyJ > W <*P = max {w a p, W 0a 

As a result, we establish 

%W N N 

|X(u,v)| <— EE (iKIk Klk + Klly, IMyJ = 

ct=l 0=1 
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6WN A 



2^ l|Ua|lvJ|Va|lv a < Z \H v \\ v \\v ' Vu > v e Vo, 



=1 



where Z = 6 l¥7V/6> > 0, W = max W a( g > 0. □ 

l<a,0<N K 

From this lemma and lemma 1, it follows that bilinear form ((75)) is sym- 
metric, continuous and coercive with constants M = M + Z and B = B respec- 
tively. In addition, due to lemmas 1, 2, and 4, functionals L(v), X(u k ,v), 
and Jg(u fe ,v) are linear and continuous in v. Therefore, there exists a unique 
solution u = u fc+1 e Vq of variational problem ((781) . 

Thus, the conditions of theorem 4 are satisfied, and we obtain the next 
proposition. 

Theorem 6. Suppose that 7 e (0;7 2 ), 72 = 2B 2 /(M + D) 2 . Then the 
sequance {u fc }, obtained by iterative method J7ff| ) - {79$ , which is equivalent to 
iterative method |^6p with bilinear form |76| ), converges strongly in Vq to the ex- 
act solution u £ Vq of nonlinear penalty variational equation pll) for unilateral 



multibody contact problem, i.e. ||u — u||^ — > 0. Moreover, the convergence 

rate in norm || • || G is linear H50\) . where q = 1 — 7 (2B — jM 2 / B) j (M + Z) , 

and the maximal rate reaches as 7 = 7 = B 2 /M 2 , M» =M + D. 

Now let us show that iterative method ([78")) - ([79")) leads to domain decom- 
position. 

Due to relationships 

1 N f 

Jg(u, v) = -- / (<f Qi g - u an - up n ) v an dS, 

a=l 0£B a J S <*P 

1 N f 

a=l gB Q ^ S °' 3 

method ((75)) - ((79")) rewrites as follows: 



N N 

E««(«o +1 . v «)+ E E / ' 

^ 1 „ 1 fl ^ D 0,vfl 



- w„ n ) v an dS 



a=l /3 SB, 



= E Za ( V «) + ^E E / ( d ^- u an-^n) «andS, (83) 
a=l a=l /3<£B a S "fi 

uW=7U^ + (l- 7 )u^ a = 1,2, ...,7V, fc = 0,l,.... (84) 

Since the common quantities of the subdomains are known from the previous 
iteration, the variational equation (|83[) splits into N variational equations in 
separate subdomains ft a . Therefore, method (|83)) - (|84)) can be written in the 
following equivalent form: 

a Q (u£ +1 , v Q ) + - V / ip a pv%£v an dS = 
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= l a (v a )+g ^ / 4>ai)U k an v an dS+ 
+ \H (d a p-u k an -u k p n y v an dS, Vv a eV°, (85) 

< +1 =7fia +1 + (l-7)u^, a = 1,2, ...,7V, fc = 0, 1, ... . (86) 
Since bilinear forms a a (u Q , v a ) , JC Q (u a , v a ) = ^J^/jes Xs fj ^aSiU an v an dS 
are symmetric, continuous and coercive, and functionals i a (v a ), X Q (u*,v Q ), 
\ J2/3eB a Is g i^ a P ~ u an ~ u pn^ v an dS are linear and continuous in v a , it 
follows that there exists a unique solution u* = u„ +1 G V„ of each varia- 
tional equation (1551 . Furthermore, it is obvious to see that the unique solu- 
tion u* of variational equation (1751 . i.e. (1551) . takes the form u* = u fc+1 = 
(u*, uj-j, Ujy) G Vq. Therefore, the solution of variational equation (|75|) is 
equivalent to the solution of iV variational equations (|85p in separate subdo- 
mains, and iterative processes (|78l) - (l79l) and ((85)) - ([86]) are equivalent. 
Now, consider iterative method (155]) - (|55|) in more detail. 
In each iteration of this method we have to solve in parallel N variational 
equations (1851) . which correspond to some elasticity problems in subdomains 
with prescribed Robin boundary conditions on possible contact areas: 

5 ^8» + ^«an7 fl = i d ^-<n-4nY / e + ^U k an /6 OnS a p. (87) 

Here 0^+1 are unknown normal stresses on possible contact areas S a p. 

Therefore, iterative method (|85p - ((55)) refers to parallel Robin-Robin type 
domain decomposition schemes. 

Since domain decomposition method (155)) - ([561 and iterative method (1751) - 
(175)) are equivalent, the convergence theorem 6 also holds for the method 
055]) - (EHD- 

Taking different characteristic functions ip a p in (|85p . i.e. different subareas 
S^p of possible contact zones S a p, we can obtain different particular cases of 
domain decomposition method ((55)) - (f55")) . 

Thus, taking ^ a| g(x) = 1, Va, /3, i.e. S 1 ^ = S^, we get domain decomposi- 
tion scheme with Robin boundary conditions on whole possible contact areas: 

+ U L a + n/ = ( rf <*/3 ~ U k an - U k 0n y /0 + U k an /9 On S aP . 

Therefore, we have named this domain decomposition method as full parallel 
Robin-Robin domain decomposition scheme [T5J [TH] ■ 

Taking ip a p(x) = 0, Va, /3, i.e. S^r = 0, we get parallel Neumann-Neumann 
domain decomposition scheme [THl 13 [T5J QI5] : 

a a (u k+1 7 v Q ) = Z Q (v Q ) + i / (d^-u^-u^) ^,^5, Vv a G V„ , 

(88) 

u tt +1 =7u Q +1 + (l- 7 )u Q , a =1,2,..., TV, fc = 0,l,.... (89) 
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In each step fc of this scheme we have to solve in parallel N variational 
equations (I88p . which correspond to elasticity problems in sub-domains with 
Neumann boundary conditions on possible contact areas: 

° k c$n = °t{in = (<k/3 -u k an -u k 0n )~ / 6 on S af3 . 

Note, that in the most general case, we can choose functions ip a /3, i-e. sur- 
faces S^p, differently for each a, (3. 

Moreover, we can choose functions ijj a p differently at each iteration k, i.e. 

^(x)=^(x) = {0, xeS Asyv{L xesy, (90) 

where S k p C S a p, k = 0, 1, ... are some given subareas of possible contact zones 
S a p, a = 1, 2, N, f3 € B a . 

As a result we obtain following nonstationary Robin-Robin type domain 
decomposition scheme 




= Uv Q ) + i Yl [ ( d ^-u k an -u k p n ) v an dS, Vv a eV°, (91) 

u^^u^ + fl-TX, a = 1,2, ...,7V, fc = 0,l,.... (92) 

This domain decomposition scheme is equivalent to nonstationary iterative 
method (|7T|) with bilinear forms 

G fe (u,v) = A(u,v) +X fe (u,v), u,ve^ , fc = 0,l,..., (93) 

where 

X k (u,v) = ^ I (u an v an ip k p + Ui3 n V0 n <il)p a ) dS, u,v G V , k = 0, 1, ... . 

(94) 

Bilinear forms X k (u, v) are symmetric, nonnegative, and continuous with 
constant Z > 0. Hence, bilinear forms (1§3"|) satisfy properties (jT2"j) - (f?4")) . where 
M = M + Z, B = B. 

It is obvious to see that condition (1751) m theorem 5 for bilinear forms (l9"3|) 
is equivalent to the following condition 

(3fc € {0, 1, ...» (Vfc > ko) {X k (u, u) > X k+1 (u, u)} , 

which by-turn is equivalent to the condition 

(3fc e {0, 1, ...}) (Vfc > ko) (Va) (V/3 G B a ) (Vx € S a/3 ) {^(x) > ^ x (x)} . 

(95) 

Therefore, from theorem 5 we obtain the next proposition about conver- 
gence of nonstationary domain decomposition scheme (|9H - (|92l) . 

Theorem 7. Suppose that 7 £ (0572), 72 = 2B 2 /(M + D) 2 , and functions 
V'aft satisfy property \95\l. Then the sequence {u fc }, obtained by nonstationary 
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domain decomposition scheme Wl\) - A9ty) , converges strongly in Vo to the exact 
solution of nonlinear penalty variational equation |^7]) for unilateral multibody 
contact problem, i.e. u fc — u ,, — > 0. 

Now let us consider a particular case of domain decomposition method (|9Tj) - 
([9"2"]l . In each iteration fc let us choose functions tjr^g as follows [TO1 [T51 [TO) 129]: 

^(x) - Xa/J (x) - | ^ ^ (x) _ uL(x) _ u , Jx , } < Q , (96) 

where x G S^, x' = P(x) G • Then, taking into consideration that 
(d aP - u k an - vf^A = (d a p - u k an - u k pr ^j x k aj3 , we obtain the method PH 

a a (u k+1 , v Q ) + i V / xlp (w^ 1 ~ - u k Pn )) v an dS = l a (v a ), (97) 

< +1 =7u" +1 + (l-7)ui a = l,2,...,JV, fc = 0,l,.... (98) 

In each step k of this method we have to solve in parallel N variational 
equations (|97|) . which correspond to elasticity problems in subdomains with 
prescribed displacements d a /3 — u k n through penalty on some subareas of pos- 
sible contact zones S a p. Therefore, we can conventionally name this method as 
nonstationary parallel Dirichlet-Dirichlet domain decomposition scheme. 

Note, that in some particular cases of penalty Robin-Robin domain decom- 
position method (|9Tj) - ((92)) (for example ip k p{x.) = 1, Va, (3 or ^^(x) = x^(x), 
Va, f3) we can pass the limit to 9 — > 0, and obtain domain decomposition schemes 
without penalty [T9l I30j . 

The advantages of proposed domain decomposition schemes are their sim- 
plicity, and the regularization of the original contact problem because of the 
use of penalty term. These domain decomposition schemes have only one itera- 
tion loop, which deals with domain decomposition and nonlinearity of unilateral 
contact conditions. 

The disadvantage of proposed methods is that we have to choose a penalty 
parameter. 

Finally, let us say that iterative methods (jH)]) and (|71D for the solution 
of nonlinear variational equations are rather general. From these methods, 
besides parallel Robin-Robin type domain decomposition schemes ([55)1 - ([55]) 
and (jnj) - ([52"]) , we also can obtain other different particular iteratve methods 
for the solution of penalty variational equation of unilateral multibody contact 
problems, which do not lead to domain decomposition. 

Thus, taking bilinear form G (u, v) = A (u, v) + X (u, v), 
X ( U > V ) = I E{a, 0} e Q f Sa p ( u « n+u Pn ){v an + v f}n ) dS , 
and iterative parameter 7 = 1 in (|46[) . we obtain iterative method for the solu- 
tion of multibody unilateral contact problems, which can be viewed as a gener- 
alization of penalty iteration method, proposed in [27] for the solution of crack 
problems with nonpenetration conditions. 

Taking in method (|7T| bilinear forms G k in each step k as follows 

G fe (u,v) = A(u,v)+X k (u,v), u,vgV , 
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where 

X fe (u,v) = i ^ I {u a n + Upn)(Van + V0 n ) X%dS, U, V € Vq, 
{a, 0}eQ Js «(> 

X ^ W - \ 1, rf^(x) - u*„(x) - u* n (x') < ' X G ^' X - ^ W G i/3Q ' 

and iterative parameter 7=1, we obtain iterative method, which can be viewed 
as active set method, i.e. semi-smooth Newton method, for unilateral multibody 
contact problems. For one-body contact problems active set method was intro- 
duced in |20j . The convergence theorem for active set method, as a variant of 
semi-smooth Newton method, was proved in |31j . 

Note, that bilinear forms X (u, v) and X k (u 7 v) are symmetric, nonnegative 
and coercive. 

However, these two particular cases of methods (l4"6l and (|7ip do not lead to 
the domain decomposition. 

In the next section we investigate numerical efficiency of proposed by us 
penalty parallel Robin-Robin type domain decomposition schemes (|85j) - ((86)) 
and © " &■ 



7 Numerical investigations 

We provide numerical investigation of proposed domain decomposition schemes 
for plane unilateral contact problems of two elastic bodies ft a C R 2 , a = 1,2. 
For the numerical solution of linear variational problems in subdomains, we use 
finite element method (FEM) with linear and quadratic triangular elements. 

At first let us compare the convergence rates of different particular domain 
decomposition schemes. 

Consider contact problem of two transversally isotropic elastic bodies £l a , 
a = 1, 2 with the plane of isotropy, parallel to the plane X2 = (fig. 2). 

The material properties of the bodies are: E a /E' a — 2, G a /G' a =2, v a = 
v' a = 0.3, where E a , v a , and G a - are the elasticity modulus, Poisson's ratio, 
and the shear modulus for the body VL a in the plane of isotropy, and E' a , v' a , 
G' a are these constants in the orthogonal direction, a = 1, 2. 

The length and height of each body is the same, and is equal to 4 b. The 
distance between bodies before the deformation is c?i2(x) = rx\/b 2 , the com- 
pression of the bodies is A w 2.154434 r, r = 10~ 3 b, and the possible contact 

area is S12 = jx = (xi, x 2 ) T : Xx € [0; 2 6], X2 = 4&j. 

The problem is solved by parallel Robin-Robin domain decomposition schemes, 
using FEM with 1595 quadratic triangular elements in each body (15 elements 
on each side of possible contact area S12). 

We take the penalty parameter in the form 9 = 4:bc(l/E[ + 1/E 2 ), c = 0.05, 
where c is dimensionless penalty coefficient. The initial guesses for displacements 
ti°„, a — 1,2 are chosen using the bar model [T9] : 

. . 4 6c \d12M - A]~ „ . , 46c[rf i2 (x) - A]~ . 
<(x) = ^ +c) , x e S 12) ,L(x) = 0E' 2 ll + c) +A ' X G 521 ' 

We use the following stopping criterion for domain decomposition schemes 
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Fig. 2. Unilateral contact of two transversaly-isotropic bodies 




h k c£ -ulX/\K + n\\ 2 <Zu, a = l,2,...,N, (99) 

where || Wq, T i 1 1 2 = \ZSj Kntx 1 )] 2 is the discrete norm, x- 7 6 6*12 are the finite 
element nodes on the possible contact area, and e u > is the relative accuracy 
for displacements. 

Fig. 3 shows the approximations to dimensionless normal contact stress 
o-*(^i,x 2 ) = (T 12 n(xi,x 2 )/ |<7i2n(0, x 2 )\, x 2 = 4 6, (xi^^ e Sv2, obtained 
by parallel Neumann-Neumann scheme ([55]) - (1551 (tp\ 2 (x) = t/i 2 i (x) = 0) at it- 
erations k — 1, 2, 4, 21 (curves 1-4) for optimal iterative parameter 7 = 0.173 
and accuracy e u = 1CP 3 . The dashed curve represents the exact solution for 
two half-spaces, obtained in |32| , Hence, the real contact area is S% 2 ~ [0j b}, 

where [y; 2] = jx = (a;i,a;2) T : Xi £ [y; z], x 2 = 4&j. 

At fig. 4 and fig. 5 the convergence rates of different particular domain 
decomposition schemes are compared. 

The total iteration number m dependence on iteration parameter 7 for ac- 
curacy e u — 10~ 3 is shown at fig. 4, and its dependence on logarithmic accuracy 
lg e u for optimal iteration parameter 7 = 7 is shown at fig. 5. 

The first curve at these figures represents the parallel Neumann-Neumann 
scheme (S\ 2 = S\\ — 0, ^12 (x) = V>2i( x ) = 0), curves 2, 3, 4 and 5 correspond to 
the parallel Robin-Robin schemes ([55)1 - (|55| with S\ 2 — S 21 equal to [0; 0.5], 
[0; 1], [0; 1.5], and [0; 2] = = S 12 , ^ 12 (x) = V21 (x) = 1) respectively. 
Curve 3 also represents the nonstationary parallel Dirichlet-Dirichlet scheme 

(SZ) -(23). 

The optimal iterative parameters 7 for the schemes represented at curves 1— 
5 are 7 = 0.173, 0.39, 0.72, 0.85 and 0.92 respectively. For 7 = 7 and accuracy 
e u = 10~ 3 these schemes converge in 21, 11, 5, 11 and 14 iterations. 
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Fig. 3. Dimensionless normal contact stress at different iterations 
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Thus, the convergence rate of stationary Robin-Robin domain decomposi- 
tion schemes is linear. The parallel Robin-Robin scheme (|55j) - (|86p with the 
surfaces Sl 2 , S 21 most closed to the real contact area (Sl 2 = S 21 ~ 5 12 ~ 
[0; &]), and the nonstationary parallel Dirichlet-Dirichlet scheme (|97|) - ([55)1 
(ipi2 — ^21 = X12)' which are represented at curve 3, have the highest con- 
vergence rates. These two schemes also have the widest range from which the 
iterative parameter 7 can be chosen. The convergence rate of parallel Neumann- 
Neumann scheme (Sl 2 = S 21 = 0), which is represented at curve 1, is the most 
slow. 

Now let us investigate the convergence of penalty method and its dependence 
on finite element discretization. 

Consider the unilateral contact problem of two isotropic bodies fij and fl 2 , 
one of which has a groove (fig. 6). 

The bodies are uniformly loaded by normal stress with intensity q. Each 
body has length I and height h, and the grove has length b. 

The material properties of the bodies are the same: E\ = E 2 = E, v\ = 
v 2 = v = 0.3. The distance between bodies before the deformation is dia(x) = 

r{ [1 — (xi — l) 2 /b 2 ] + Y^ 2 , where r = 0.056, y + = max{0,y}. The possible 
contact area is S\ 2 = |x = (a;i,a;2) T : Xi £ [0; I], x 2 — hX. 

The problem is solved by nonstationary parallel Dirichlet-Dirichlet domain 
decomposition scheme (|97p - (|98[) with finite element approximations on trian- 
gles. 

The penalty parameter is taken as follows 

2 

6 = ch (1 - v a f J E a , (100) 

Q = l 

where c is dimensionless penalty coefficient. Initial guesses for displacements 
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Fig. 6. Unilateral contact of two bodies with a groove 



are chosen using the bar model [19 : 



t (x) 



«5„(x) 



= eft (l-i/g) [d 12 (x)-A] 
6»Si (1 + c) 

eft [d 12 (x)-A]~ 
9E 2 (1 + c) 



xeS, 



A, xe& 



where A = q 9 (1 + c) /c. We use (|99p as stopping criterion for iterative process. 

For iterative parameter 7 6 [0.45; 0.65], accuracy e„ = 10~ 3 , penalty coeffi- 
cients c and finite element meshes considered below, parallel Dirichlet-Dirichlet 
scheme, applied to solve this problem, converges in 2-15 iterations. 

We investigate the dependence of the quality of numerical solution, obtained 
by this scheme, on penalty parameter and finite element mesh. 

Plots at fig. 7 represent approximations to dimensionless contact stress <r* (x) = 
<?i2 n( x ) /E, x 6 Si 2 for bodies with size h — I = 8 b and external load q = 0.01 E, 
obtained by Dirichlet-Dirichlet scheme for different dimensionless penalty coef- 
ficients c at fixed finite element mesh with 64 linear triangular finite elements 
on each side of possible contact area S12. Curves 1-4 correspond to c = 0.1, 
0.05, 0.01 and 0.0025 respectively. 

Plots at fig. 8 represent approximations to ff*(x), x G S12 for bodies with 
length I — 8 6, height h — 2b, and external load q — 0.0075 E, obtained for 
differrent dimensionless penalty coefficients c and different finite element meshes. 

Curves 1 and 2 at this figure correspond to er* for dimensionless penalty 
coefficients c = 0.1 and c = 0.01 respectively at finite element mesh with 32 
linear triangular elements on each side of possible unilateral contact area S±2- 
Curves 3 and 4 correspond to a* for c = 0.1 and c = 0.01 respectively, but 
for finite element mesh with 64 linear triangular elements on each side of Si2- 
Dashed curve at this figure and at fig. 7 represents the exact solution, obtained 
in [33 for the contact of two half-spaces. Here we see that in spite of the 
solution, obtained for penalty coefficient c = 0.1 at mesh with 32 finite elements 
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Fig. 7. Normal contact stress <r* for different penalty coefficients (at fixed 

mesh) 
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on each side of £12 (curve 1 at fig. 8), the solution obtained for lower penalty 
coefficient c = 0.01 at the same mesh become instable. But if we refine finite 
element mesh twice for the penalty coefficient c = 0.01, then the influence of 
errors on the perturbation of initial data will decrease, and we shall obtain much 
better approximation to the exact solution (curve 4 at fig. 8). 

Hence, we conclude that for obtaining a nice approximation to the solu- 
tion, we need to decrease penalty parameter and to refine finite element mesh 
simultaneously. 

8 Conclusions 

For the solution of unilateral multibody contact problems of elasticity we have 
proposed on continuous level a class of parallel Robin-Robin type domain de- 
composition schemes, which are based on penalty method for variational in- 
equalities and some stationary or nonstationary iterative methods for nonlinear 
variational equations. In each step of these schemes we have to solve in parallel 
linear variational equations in subdomains, which correspond to some elasticity 
problems with Robin boundary conditions on possible contact areas. 

We have proved the theorems about the strong convergence and stability 
of these schemes, and we have shown that the convergence rate of stationary 
Robin-Robin schemes in some energy norm is linear. 

Numerical analysis of proposed domain decomposition schemes has been 
made for plane two-body contact problems using linear and quadratic finite 
element approximations on triangles. Convergence rates of different particular 
domain decomposition schemes have been compared and their dependence on 
iterative parameter 7 has been investigated. The penalty parameter and mesh 
refinement influence on the numerical solution have been examined. Numerical 
experiments have confirmed the theoretical results about convergence of these 
domain decomposition schemes. 

The advantages of proposed domain decomposition schemes are their sim- 
plicity, and the regularization of the original contact problem because of the use 
of penalty term. These domain decomposition schemes have only one iteration 
loop, which deals with domain decomposition and nonlinearity of unilateral con- 
tact conditions. In addition, these methods are obtained on continuous level, 
and they do not depend on any discretization technique. 

The generalization of proposed domain decomposition schemes for contact 
problems of nonlinear elastic bodies are possible. 

Moreover, we have built similar domain decomposition schemes for ideal 
mechanical contact problems (see our work [T"9"]). 

The disadvantage of proposed methods is that we have to choose a penalty 
parameter. 
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